function [ssT_draw, rho_inf] = regime_prob(para,data)

global mprior psel 
 
if mprior == 1 
 % calculates log-likelihood based on Kalman filter  

alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a  = para(8);
rho_f  = para(9);
rho_i  = para(10);
sd_a   = para(11);
sd_f   = para(12);
sd_i   = para(13);
piess  = para(14);
gy     = para(15);
lnA_0  = para(16);
yst    = para(17);
p11    = para(18);
p22    = para(19);

                 
iss    = gy-log(beta1) + piess;

        
%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------

A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];


eu = 0; 



A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;kappa1;kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];
 y1a = solvec_a(1);
 y2a = solvec_a(2);
 pi1a = solvec_a(3);
 pi2a = solvec_a(4);
 i1a = solvec_a(5);
 i2a = solvec_a(6);
 y1f = solvec_f(1);
 y2f = solvec_f(2);
 pi1f = solvec_f(3);
 pi2f = solvec_f(4);
 i1f = solvec_f(5);
 i2f = solvec_f(6);
 y1i = solvec_i(1);
 y2i = solvec_i(2);
 pi1i = solvec_i(3);
 pi2i = solvec_i(4);
 i1i=solvec_i(5);
 i2i=solvec_i(6);

%----------------------------------------------------------------------
% compute term structure
%----------------------------------------------------------------------
rho      = diag([rho_a;rho_f;rho_i]);
sd       = [sd_a;sd_f;sd_i];
sd1 = sd; sd2 = sd; 
sd_i1=sd_i; sd_i2=sd_i; 


B1 = zeros(1,3);
B2 = zeros(1,3);


B1(1,:) =  [i1a, i1f, i1i];
B2(1,:) =  [i2a, i2f, i2i];



% Augment the state vector 
rho = diag([diag(rho);1]);
rho(4,1) = rho_a; 
TC  = [0;0;0;gy]; 
sd1  = [sd1;sd_a];
sd2  = [sd2;sd_a]; 
Vs1 = diag([sd_a^2/(1-rho_a^2);sd_f^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0]);
Vs2 = diag([sd_a^2/(1-rho_a^2);sd_f^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);

% Inflation Persistence 
sumw1      = pi1a^2*Vs1(1,1)+pi1f^2*Vs1(2,2)+pi1i^2*Vs1(3,3); 
sumw2      = pi2a^2*Vs1(1,1)+pi2f^2*Vs1(2,2)+pi2i^2*Vs1(3,3); 
weight11   = (pi1a^2)*Vs1(1,1)/sumw1; 
weight12   = (pi1f^2)*Vs1(2,2)/sumw1; 
weight21   = (pi2a^2)*Vs1(1,1)/sumw2; 
weight22   = (pi2f^2)*Vs1(2,2)/sumw2; 
rho_inf1   = rho_a*weight11+rho_f*weight12+rho_i*(1-weight11-weight12);  
rho_inf2   = rho_a*weight21+rho_f*weight22+rho_i*(1-weight21-weight22); 
rho_inf    = [rho_inf1;rho_inf2]; 

% system matrices for the macro part 
mac1 = [y1a y1f y1i 1; pi1a pi1f pi1i 0]; 
mac2 = [y2a y2f y2i 1; pi2a pi2f pi2i 0]; 

% Markov Chain transition probability 
trans = [ p11, 1-p11;1-p22, p22];  
trans = trans';
probs = [(1-p22)/(2-p11-p22);(1-p11)/(2-p11-p22)];

Vs  = probs(1)*Vs1+probs(2)*Vs2; 
% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 
nhist = 3;              
nxmax = 2^nhist; 

Az      = [ yst ; piess ; iss]; Az1 = Az; Az2 = Az; 
Bz      = [ mac1 ; B1(1,:) 0]; 
Bz1     = [ mac1 ; B1(1,:) 0];    
Bz2     = [ mac2 ; B2(1,:) 0]; 

% The matrix for storing information from the forward simulation 
probsT = zeros(nobs,2); 
ssT_draw = zeros(nobs,2); 
probsT(1,:) = probs';


% condition on the first observation 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 

llik_c   = zeros(nobs,1);       % log likelihood contribution
nu        = data(1,:)'-(Az+Bz*(TC+rho*s00)); 
Ft        = Bz*Vs*Bz'; 
llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu'*inv(Ft)*nu;

s00      = s00+(Vs*Bz')*inv(Ft)*nu; 
ltrvecP1 = (ltrvec(Vs-Vs*Bz'*inv(Ft)*(Vs*Bz')'))';
probx    = 1;
s1t = s00'; 
t=2;
% log likelihood by Kalman filter
while t<nobs+1
nx    = size(s1t,1);   
probs = trans*probs; 
probx2 = kron(probs,probx); 

llik_ht  = zeros(2*nx,1);      
s0t       = zeros(2*nx,ns);    
ltrvecP0 = zeros(2*nx,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 

i = 1; 
 while i<=nx; 

     
   % forecasting
   
   alphahat = TC+rho*s1t(i,:)';
   P1t      = ltrmat((ltrvecP1(i,:))');
   Phat1    = rho*P1t*rho'+(diag(sd1))^2;
   Phat2    = rho*P1t*rho'+(diag(sd2))^2; 
   yhat1    = Az1 + Bz1*alphahat;
   yhat2    = Az2 + Bz2*alphahat; 
   
   nu1      = data(t,:)-yhat1'; 
   nu2      = data(t,:)-yhat2'; 
   
   Ft1      = Bz1*Phat1*Bz1';
   Ft2      = Bz2*Phat2*Bz2'; 
   
   yhat     = yhat + probx2(i)*yhat1+probx2(nx+i)*yhat2;  
   Ft       = Ft+probx2(i)*(Ft1+yhat1*yhat1')+probx2(nx+i)*(Ft2+yhat2*yhat2'); 
   
   % compute likelihood 
   
   llik_ht(i,1)    = -0.5*nz*log(2*pi)-0.5*log(det(Ft1))-0.5*nu1*inv(Ft1)*nu1'; 
   llik_ht(nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft2))-0.5*nu2*inv(Ft2)*nu2'; 
   
   % updating step 
   s0t(i,:)    = (alphahat+Phat1*Bz1'*inv(Ft1)*nu1')';
   s0t(nx+i,:) = (alphahat+Phat2*Bz2'*inv(Ft2)*nu2')'; 
   
   ltrvecP0(i,:)=(ltrvec(Phat1-Phat1*Bz1'*inv(Ft1)*Bz1*Phat1'))'; 
   ltrvecP0(nx+i,:) = (ltrvec(Phat2-Phat2*Bz2'*inv(Ft2)*Bz2*Phat2'))'; 
   
   i=i+1;
   
 end 
   
   Ft = Ft -yhat*yhat'; 
   nuhat(t,:) = data(t,:)-yhat';
   
   % llik_c computation 
   
   llik_htmax = max(llik_ht); 
   llik_c(t) = llik_c(t) + llik_htmax + log(probx2'*exp(llik_ht-llik_htmax)); 
   
   % updating mixture probabilities and the state probabilities 
   
   probx2 = (probx2.*exp(llik_ht-llik_htmax))/(probx2'*exp(llik_ht-llik_htmax));
   probs  = [sum(probx2(1:nx));sum(probx2(nx+1:2*nx))]; 
   
   % collapse the state vector 
   
   zeroprob = (probx2<=1e-5); 
   
   if sum(zeroprob) >= (2*nx-nxmax)
   % simply delete the zero probability states 
   s1t = s0t(zeroprob==0,:);
   ltrvecP1 = ltrvecP0(zeroprob==0,:); 
   probx = probx2(zeroprob==0); 
   
   else

       % collapse some of the states
   selmat = kron(eye(nx),[1 1]); 
   sumprobx2 = selmat*probx2;
   sumzeroprob = selmat*zeroprob;
   probscore  = sumprobx2+sumzeroprob;
   [B,index] =sortrows(probscore);
   collind   = index(1:2*nx-nxmax-sum(zeroprob)); 
   
   
   s1t = zeros(nxmax,ns);
   ltrvecP1 = zeros(nxmax,ntr); 
   probx = zeros(nxmax,1); 
   j=1; 
   
   i=1;
   while i<=nx; 
   
       i2 = (i-1)*2+1; 
       
       if sumzeroprob(i)>0; 
        % no merge, just eliminate 
        
         if zeroprob(i2)==0 
            s1t(j,:)=s0t(i2,:);
            ltrvecP1(j,:)=ltrvecP0(i2,:);
            probx(j,:) = probx2(i2); 
            j=j+1;
         end 
         
         if zeroprob(i2+1)==0
            s1t(j,:)=s0t(i2+1,:);
            ltrvecP1(j,:)=ltrvecP0(i2+1,:); 
            probx(j,:)=probx2(i2+1);
            j=j+1;
         end 
         
       else
           
         x=i==collind;
         x=sum(x);
         if x==1
         % collapse density 
            probx(j) = probx2(i2)+probx2(i2+1); 
            s1t(j,:) = (probx2(i2)/probx(j))*s0t(i2,:)+(probx2(i2+1)/probx(j))*s0t(i2+1,:);
            ltrvecP1(j,:)=ltrvec((probx2(i2))/(probx(j))*(ltrmat(ltrvecP0(i2,:)')+s0t(i2,:)'*s0t(i2,:)) + (probx2(i2+1))/(probx(j))*(ltrmat(ltrvecP0(i2+1,:)')+s0t(i2+1,:)'*s0t(i2+1,:))-s1t(j,:)'*s1t(j,:))';
            j=j+1;
         else
          % don't collapse density 
             s1t(j,:)=s0t(i2,:);
             ltrvecP1(j,:) = ltrvecP0(i2,:); 
             probx(j,:) = probx2(i2); 
             j=j+1; 
             s1t(j,:)=s0t(i2+1,:);
             ltrvecP1(j,:) = ltrvecP0(i2+1,:);
             probx(j,:) = probx2(i2+1);
              j=j+1;
         end

       end
       i=i+1;
   end
   end
   probx = probx./sum(probx);
   nnx = size(s1t,1);
      if nnx>nxmax;
       disp('error : nx is greater than nxmax')
       t=nobs+1; 
      end
      % store information for future use 
   probsT(t,:) = probs'; 
t=t+1;
end

logl = sum(llik_c);

% Backward Simulation 
t = nobs;
probs1   = probsT(t,1);
ssT_draw(t,:) = [probs1, 1-probs1];

while t>1 
 t = t-1; 
 probs1 = ssT_draw(t+1,:)*trans(:,1)*probsT(t,1)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 ssT_draw(t,:) = [probs1, 1-probs1]; 
end   

elseif mprior ==2  
 alpha1 = para(1);
gamma1 = para(2);
beta1  = para(3);
kappa1 = para(4);
tau1   = para(5);
rho_a1 = para(6);
rho_f  = para(7);
rho_i  = para(8);
sd_a1   = para(9);
sd_f1   = para(10);
sd_i1   = para(11);
sd_a2   = para(12);
sd_f2   = para(13);
sd_i2   = para(14);
piess  = para(15);
gy     = para(16);
lnA_0  = para(17);
yst    = para(18);
p11    = para(19);
p22    = para(20);


bound  = zeros(length(para),2);
bound(17:18,1) = -100; 
bound(:,2) = [100;100;0.9999;100;100;0.9999;0.9999;0.9999;100;100;100;100;100;100;100;100;100;100;0.9999;0.9999]; 



crit1  = para-bound(:,1); crit2=bound(:,2)-para; 
                 if min(crit1)<0 || min(crit2)<0
logl = 1000000; 
                 return 
                 else
iss    = gy-log(beta1) + piess;

%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------

alpha2 = alpha1; gamma2=gamma1; rho_a2 = rho_a1; 

A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];

eu = 0; 
if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  DL determinacy condition violated  *************')
    eu = 1;
end


                   if eu==0
A_Aa_mat = [1-p11*rho_a1, -(1-p11)*rho_a2, -(p11*rho_a1)/tau1, -((1-p11)*rho_a2)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a1, 1-p22*rho_a2, -((1-p22)*rho_a1)/tau1, -(p22*rho_a2)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a1, -beta1*(1-p11)*rho_a2, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a1, 1-beta1*p22*rho_a2, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
solvec_a = inv(A_Aa_mat) * [(p11*rho_a1+(1-p11)*rho_a2)/tau1;((1-p22)*rho_a1+p22*rho_a2)/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;kappa1;kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];
 y1a = solvec_a(1);
 y2a = solvec_a(2);
 pi1a = solvec_a(3);
 pi2a = solvec_a(4);
 i1a = solvec_a(5);
 i2a = solvec_a(6);
 y1f = solvec_f(1);
 y2f = solvec_f(2);
 pi1f = solvec_f(3);
 pi2f = solvec_f(4);
 i1f = solvec_f(5);
 i2f = solvec_f(6);
 y1i = solvec_i(1);
 y2i = solvec_i(2);
 pi1i = solvec_i(3);
 pi2i = solvec_i(4);
 i1i=solvec_i(5);
 i2i=solvec_i(6);

%----------------------------------------------------------------------
% compute term structure
%----------------------------------------------------------------------
rho1      = diag([rho_a1;rho_f;rho_i]);
rho2      = diag([rho_a2;rho_f;rho_i]);
lambda_a1 = tau1*y1a+1+pi1a;
lambda_a2 = tau1*y2a+1+pi2a;
lambda_f1 = tau1*y1f+pi1f;
lambda_f2 = tau1*y2f+pi2f;
lambda_i1 = tau1*y1i+pi1i;
lambda_i2 = tau1*y2i+pi2i;
lambda1   = [lambda_a1;lambda_f1;lambda_i1];
lambda2   = [lambda_a2;lambda_f2;lambda_i2];

sd1 = [sd_a1;sd_f1;sd_i1]; 
sd2 = [sd_a2;sd_f2;sd_i2];  
 



B1 = zeros(1,3);
B2 = zeros(1,3);

B1(1,:) =  [i1a, i1f, i1i];
B2(1,:) =  [i2a, i2f, i2i];

% Augment the state vector 
rho1 = diag([diag(rho1);1]);
rho2 = diag([diag(rho2);1]); 
rho1(4,1) = rho_a1;
rho2(4,1) = rho_a2; 
TC  = [0;0;0;gy]; 
sd1  = [sd1;sd_a1];
sd2  = [sd2;sd_a2]; 
Vs1      = diag([sd_a1^2/(1-rho_a1^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0]);
Vs2      = diag([sd_a2^2/(1-rho_a2^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);
% Inflation Persistence 
sumw1      = pi1a^2*Vs1(1,1)+pi1f^2*Vs1(2,2)+pi1i^2*Vs1(3,3); 
sumw2      = pi2a^2*Vs2(1,1)+pi2f^2*Vs2(2,2)+pi2i^2*Vs2(3,3); 
weight11   = (pi1a^2)*Vs1(1,1)/sumw1; 
weight12   = (pi1f^2)*Vs1(2,2)/sumw1; 
weight21   = (pi2a^2)*Vs2(1,1)/sumw2; 
weight22   = (pi2f^2)*Vs2(2,2)/sumw2; 
rho_inf1   = rho_a1*weight11+rho_f*weight12+rho_i*(1-weight11-weight12);  
rho_inf2   = rho_a2*weight21+rho_f*weight22+rho_i*(1-weight21-weight22); 
rho_inf    = [rho_inf1;rho_inf2]; 


% system matrices for the macro part 
mac1 = [y1a y1f y1i 1; pi1a pi1f pi1i 0]; 
mac2 = [y2a y2f y2i 1; pi2a pi2f pi2i 0]; 

% Markov Chain transition probability 
trans = [ p11, 1-p11;1-p22, p22];  
trans = trans'; 
probs = [(1-p22)/(2-p11-p22);(1-p11)/(2-p11-p22)];


Vs1 = probs(1)*diag([sd_a1^2/(1-rho_a1^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0])+probs(2)*diag([sd_a2^2/(1-rho_a2^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);

% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho1,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 
nhist = 3;              
nxmax = 2^nhist; 

Az      = [ yst ; piess ; iss]; Az1 = Az; Az2 = Az; 
Bz      = [ mac1 ; B1(1,:) 0]; 
Bz1     = [ mac1 ; B1(1,:) 0];    
Bz2     = [ mac2 ; B2(1,:) 0]; 


% The matrix for storing information from the forward simulation 
probsT = zeros(nobs,2); 
ssT_draw = zeros(nobs,2); 
probsT(1,:) = probs';


% condition on the first observation 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 

llik_c   = zeros(nobs,1);       % log likelihood contribution
nu        = data(1,:)'-(Az+Bz*(TC+rho1*s00)); 
Ft        = Bz*Vs1*Bz'; 
llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu'*inv(Ft)*nu;

s00      = s00+(Vs1*Bz')*inv(Ft)*nu; 
ltrvecP1 = (ltrvec(Vs1-Vs1*Bz'*inv(Ft)*(Vs1*Bz')'))';
probx    = 1;
s1t = s00'; 
t=2;
% log likelihood by Kalman filter
while t<nobs+1
nx    = size(s1t,1);   
probs = trans*probs; 
probx2 = kron(probs,probx); 

llik_ht  = zeros(2*nx,1);      
s0t       = zeros(2*nx,ns);    
ltrvecP0 = zeros(2*nx,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 

i = 1; 
 while i<=nx; 

     
   % forecasting
   
   alphahat1 = TC+rho1*s1t(i,:)';
   alphahat2 = TC+rho2*s1t(i,:)';
   P1t      = ltrmat((ltrvecP1(i,:))');
   Phat1    = rho1*P1t*rho1'+(diag(sd1))^2;
   Phat2    = rho2*P1t*rho2'+(diag(sd2))^2;
   yhat1    = Az1 + Bz1*alphahat1;
   yhat2    = Az2 + Bz2*alphahat2; 
   
   nu1      = data(t,:)-yhat1'; 
   nu2      = data(t,:)-yhat2'; 
   
   Ft1      = Bz1*Phat1*Bz1';
   Ft2      = Bz2*Phat2*Bz2'; 
   
   yhat     = yhat + probx2(i)*yhat1+probx2(nx+i)*yhat2;  
   Ft       = Ft+probx2(i)*(Ft1+yhat1*yhat1')+probx2(nx+i)*(Ft2+yhat2*yhat2'); 
   
   % compute likelihood 
   
   llik_ht(i,1)    = -0.5*nz*log(2*pi)-0.5*log(det(Ft1))-0.5*nu1*inv(Ft1)*nu1'; 
   llik_ht(nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft2))-0.5*nu2*inv(Ft2)*nu2'; 
   
   % updating step 
   s0t(i,:)    = (alphahat1+Phat1*Bz1'*inv(Ft1)*nu1')';
   s0t(nx+i,:) = (alphahat1+Phat2*Bz2'*inv(Ft2)*nu2')'; 
   
   ltrvecP0(i,:)=(ltrvec(Phat1-Phat1*Bz1'*inv(Ft1)*Bz1*Phat1'))'; 
   ltrvecP0(nx+i,:) = (ltrvec(Phat2-Phat2*Bz2'*inv(Ft2)*Bz2*Phat2'))'; 
   
   i=i+1;
   
 end 
   
   Ft = Ft -yhat*yhat'; 
   nuhat(t,:) = data(t,:)-yhat';
   
   % llik_c computation 
   
   llik_htmax = max(llik_ht); 
   llik_c(t) = llik_c(t) + llik_htmax + log(probx2'*exp(llik_ht-llik_htmax)); 
   
   % updating mixture probabilities and the state probabilities 
   
   probx2 = (probx2.*exp(llik_ht-llik_htmax))/(probx2'*exp(llik_ht-llik_htmax));
   probs  = [sum(probx2(1:nx));sum(probx2(nx+1:2*nx))]; 
   
   % collapse the state vector 
   
   zeroprob = (probx2<=1e-5); 
   
   if sum(zeroprob) >= (2*nx-nxmax)
   % simply delete the zero probability states 
   s1t = s0t(zeroprob==0,:);
   ltrvecP1 = ltrvecP0(zeroprob==0,:); 
   probx = probx2(zeroprob==0); 
   
   else

       % collapse some of the states
   selmat = kron(eye(nx),[1 1]); 
   sumprobx2 = selmat*probx2;
   sumzeroprob = selmat*zeroprob;
   probscore  = sumprobx2+sumzeroprob;
   [B,index] =sortrows(probscore);
   collind   = index(1:2*nx-nxmax-sum(zeroprob)); 
   
   
   s1t = zeros(nxmax,ns);
   ltrvecP1 = zeros(nxmax,ntr); 
   probx = zeros(nxmax,1); 
   j=1; 
   
   i=1;
   while i<=nx; 
   
       i2 = (i-1)*2+1; 
       
       if sumzeroprob(i)>0; 
        % no merge, just eliminate 
        
         if zeroprob(i2)==0 
            s1t(j,:)=s0t(i2,:);
            ltrvecP1(j,:)=ltrvecP0(i2,:);
            probx(j,:) = probx2(i2); 
            j=j+1;
         end 
         
         if zeroprob(i2+1)==0
            s1t(j,:)=s0t(i2+1,:);
            ltrvecP1(j,:)=ltrvecP0(i2+1,:); 
            probx(j,:)=probx2(i2+1);
            j=j+1;
         end 
         
       else
           
         x=i==collind;
         x=sum(x);
         if x==1
         % collapse density 
            probx(j) = probx2(i2)+probx2(i2+1); 
            s1t(j,:) = (probx2(i2)/probx(j))*s0t(i2,:)+(probx2(i2+1)/probx(j))*s0t(i2+1,:);
            ltrvecP1(j,:)=ltrvec((probx2(i2))/(probx(j))*(ltrmat(ltrvecP0(i2,:)')+s0t(i2,:)'*s0t(i2,:)) + (probx2(i2+1))/(probx(j))*(ltrmat(ltrvecP0(i2+1,:)')+s0t(i2+1,:)'*s0t(i2+1,:))-s1t(j,:)'*s1t(j,:))';
            j=j+1;
         else
          % don't collapse density 
             s1t(j,:)=s0t(i2,:);
             ltrvecP1(j,:) = ltrvecP0(i2,:); 
             probx(j,:) = probx2(i2); 
             j=j+1; 
             s1t(j,:)=s0t(i2+1,:);
             ltrvecP1(j,:) = ltrvecP0(i2+1,:);
             probx(j,:) = probx2(i2+1);
              j=j+1;
         end

       end
       i=i+1;
   end
   end
   probx = probx./sum(probx);
   nnx = size(s1t,1);
      if nnx>nxmax;
       disp('error : nx is greater than nxmax')
       t=nobs+1; 
      end
       % store information for future use 
   probsT(t,:) = probs'; 

t=t+1;
end

logl = -sum(llik_c);
  
         % Backward Simulation 
t = nobs;
probs1   = probsT(t,1);
ssT_draw(t,:) = [probs1, 1-probs1];

while t>1 
 t = t-1; 
 probs1 = ssT_draw(t+1,:)*trans(:,1)*probsT(t,1)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 ssT_draw(t,:) = [probs1, 1-probs1]; 
end

  elseif eu>0  
logl = 1000000;
 end                  % determinacy case 
end                   %bound 

 
 
elseif mprior==4 
alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a = para(8);
rho_f = para(9);
rho_i  = para(10);
sd_a1  = para(11);
sd_f1  = para(12);
sd_i1   = para(13);
sd_a2  = para(14);
sd_f2  = para(15);
sd_i2   = para(16);
piess  = para(17);
gy     = para(18);
lnA_0  = para(19);
yst    = para(20);
p11 = para(21);
p22 = para(22);
p11_f = para(23);
p22_f = para(24);

bound  = zeros(length(para),2);
bound(19:20,1) = -100; 
bound(:,2) = [100;100;100;100;0.9999;100;100;0.9999;0.9999;0.9999;100;100;100;100;100;100;100;100;100;100;1;1;1;1]; 
     
crit1 = para-bound(:,1); 
crit2 = bound(:,2)-para; 

        if min(crit1)<0 || min(crit2)<0 
logl = -10000000; 
           return
        else


P_i = [p11 1-p11;
       1-p22 p22];
   
P_f = [p11_f 1-p11_f;
       1-p22_f p22_f];   

P = kron(P_i,P_f);

iss    = gy-log(beta1) + piess;

%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------
A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];

eu =0; 
if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  DL determinacy condition violated  *************')
    eu =1; 
end

if eu ==0 
%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------
A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;kappa1;kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];
 y1a = solvec_a(1);
 y2a = solvec_a(2);
 pi1a = solvec_a(3);
 pi2a = solvec_a(4);
 i1a = solvec_a(5);
 i2a = solvec_a(6);
 y1f = solvec_f(1);
 y2f = solvec_f(2);
 pi1f = solvec_f(3);
 pi2f = solvec_f(4);
 i1f = solvec_f(5);
 i2f = solvec_f(6);
 y1i = solvec_i(1);
 y2i = solvec_i(2);
 pi1i = solvec_i(3);
 pi2i = solvec_i(4);
 i1i=solvec_i(5);
 i2i=solvec_i(6);

%----------------------------------------------------------------------
% compute coefficients for measurement equations 
%----------------------------------------------------------------------
rho     = diag([rho_a;rho_f;rho_i]);

sd1       = [sd_a1;sd_f1;sd_i1];
sd2       = [sd_a2;sd_f2;sd_i2]; 
sd3 =sd1; sd4 =sd2; 


B1 = zeros(1,3); B2 = zeros(1,3);
B3 = zeros(1,3); B4 = zeros(1,3); 

B1(1,:) =  [i1a, i1f, i1i];
B2(1,:) =  [i1a, i1f, i1i];
B3(1,:) =  [i2a, i2f, i2i];
B4(1,:) =  [i2a, i2f, i2i];

% Augment the state vector 
rho1 = diag([diag(rho);1]);
rho2 = diag([diag(rho);1]);
rho3 = diag([diag(rho);1]);
rho4 = diag([diag(rho);1]);
rho1(4,1) = rho_a; 
rho2(4,1) = rho_a; 
rho3(4,1) = rho_a; 
rho4(4,1) = rho_a; 
TC  = [0;0;0;gy]; 
sd1  = [sd1;sd_a1]; sd2  = [sd2;sd_a2];  
sd3  = [sd3;sd_a1]; sd4  = [sd4;sd_a2]; 
Vs1 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0]);
Vs2 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);
Vs3 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0]);
Vs4 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0]);

 % Inflation Persistence 
sumw1      = pi1a^2*Vs1(1,1)+pi1f^2*Vs1(2,2)+pi1i^2*Vs1(3,3); 
sumw2      = pi1a^2*Vs2(1,1)+pi1f^2*Vs2(2,2)+pi1i^2*Vs2(3,3);
sumw3      = pi2a^2*Vs3(1,1)+pi2f^2*Vs3(2,2)+pi2i^2*Vs3(3,3); 
sumw4      = pi2a^2*Vs4(1,1)+pi2f^2*Vs4(2,2)+pi2i^2*Vs4(3,3);

weight11   = (pi1a^2)*Vs1(1,1)/sumw1; 
weight12   = (pi1f^2)*Vs1(2,2)/sumw1; 
weight21   = (pi1a^2)*Vs2(1,1)/sumw2; 
weight22   = (pi1f^2)*Vs2(2,2)/sumw2; 
weight31   = (pi2a^2)*Vs3(1,1)/sumw3; 
weight32   = (pi2f^2)*Vs3(2,2)/sumw3; 
weight41   = (pi2a^2)*Vs4(1,1)/sumw4; 
weight42   = (pi2f^2)*Vs4(2,2)/sumw4; 


rho_inf1   = rho_a*weight11+rho_f*weight12+rho_i*(1-weight11-weight12);  
rho_inf2   = rho_a*weight21+rho_f*weight22+rho_i*(1-weight21-weight22); 
rho_inf3   = rho_a*weight31+rho_f*weight32+rho_i*(1-weight31-weight32);  
rho_inf4   = rho_a*weight41+rho_f*weight42+rho_i*(1-weight41-weight42); 


rho_inf    = [rho_inf1;rho_inf2;rho_inf3;rho_inf4]; 

% system matrices for the macro part 
mac1 = [y1a y1f y1i 1; pi1a pi1f pi1i 0]; 
mac2 = [y1a y1f y1i 1; pi1a pi1f pi1i 0]; 
mac3 = [y2a y2f y2i 1; pi2a pi2f pi2i 0]; 
mac4 = [y2a y2f y2i 1; pi2a pi2f pi2i 0]; 

% Markov Chain transition probability 
trans = P';  

% compute the stationary distribution 
A = [eye(length(trans))-trans;ones(1,length(trans))]; 
probs = inv(A'*A)*(A'*[zeros(length(P),1);1]); 
Vs  = probs(1)*Vs1+probs(2)*Vs2+probs(3)*Vs3+probs(4)*Vs4; 

% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho1,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 
nhist = 2;              
nxmax = 4^nhist; 

Az      = [ yst ; piess ; iss]; 
Az1 = Az; Az2 = Az; Az3 = Az; Az4 = Az; 
Bz1     = [ mac1 ; B1(1,:) 0];    
Bz2     = [ mac2 ; B2(1,:) 0]; 
Bz3     = [ mac3 ; B3(1,:) 0];    
Bz4     = [ mac4 ; B4(1,:) 0]; 


% The matrix for storing information from the forward simulation 
probsT = zeros(nobs,4); 
ssT_draw = zeros(nobs,4); 
probsT(1,:) = probs';


% condition on the first observation : assume the first regime 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 

llik_c   = zeros(nobs,1);       % log likelihood contribution
nu        = data(1,:)'-(Az+Bz1*(TC+rho1*s00)); 
Ft        = Bz1*Vs*Bz1'; 
llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu'*inv(Ft)*nu;

s00      = s00+(Vs*Bz1')*inv(Ft)*nu; 
ltrvecP1 = (ltrvec(Vs-Vs*Bz1'*inv(Ft)*(Vs*Bz1')'))';
probx    = 1;
s1t = s00'; 
t=2;
% log likelihood by Kalman filter
while t<nobs+1
nx    = size(s1t,1);   
probs = trans*probs; 
probx2 = kron(probs,probx); 

llik_ht  = zeros(4*nx,1);      
s0t      = zeros(4*nx,ns);    
ltrvecP0 = zeros(4*nx,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 

i = 1; 
 while i<=nx; 

     
   % forecasting
   
   alphahat1= TC+rho1*s1t(i,:)';
   alphahat2= TC+rho2*s1t(i,:)';
   alphahat3= TC+rho3*s1t(i,:)';
   alphahat4= TC+rho4*s1t(i,:)';
   P1t      = ltrmat((ltrvecP1(i,:))');
   Phat1    = rho1*P1t*rho1'+(diag(sd1))^2; 
   Phat2    = rho2*P1t*rho2'+(diag(sd2))^2; 
   Phat3    = rho3*P1t*rho3'+(diag(sd3))^2; 
   Phat4    = rho4*P1t*rho4'+(diag(sd4))^2; 
   
   yhat1    = Az1 + Bz1*alphahat1;
   yhat2    = Az2 + Bz2*alphahat2; 
   yhat3    = Az3 + Bz3*alphahat3; 
   yhat4    = Az4 + Bz4*alphahat4; 
   
   
   nu1      = data(t,:)-yhat1'; 
   nu2      = data(t,:)-yhat2'; 
   nu3      = data(t,:)-yhat3'; 
   nu4      = data(t,:)-yhat4'; 
   
   Ft1      = Bz1*Phat1*Bz1';
   Ft2      = Bz2*Phat2*Bz2'; 
   Ft3      = Bz3*Phat3*Bz3'; 
   Ft4      = Bz4*Phat4*Bz4'; 
   
   yhat     = yhat + probx2(i)*yhat1+probx2(nx+i)*yhat2+probx2(2*nx+i)*yhat3+probx2(3*nx+i)*yhat4;  
   Ft       = Ft+probx2(i)*(Ft1+yhat1*yhat1')+probx2(nx+i)*(Ft2+yhat2*yhat2')+probx2(2*nx+i)*(Ft3+yhat3*yhat3')+probx2(3*nx+1)*(Ft4+yhat4*yhat4'); 
   
   % compute likelihood 
   
   llik_ht(i,1)    = -0.5*nz*log(2*pi)-0.5*log(det(Ft1))-0.5*nu1*inv(Ft1)*nu1'; 
   llik_ht(nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft2))-0.5*nu2*inv(Ft2)*nu2'; 
   llik_ht(2*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft3))-0.5*nu3*inv(Ft3)*nu3'; 
   llik_ht(3*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft4))-0.5*nu4*inv(Ft4)*nu4'; 
    
   % updating step 
   s0t(i,:)    = (alphahat1+Phat1*Bz1'*inv(Ft1)*nu1')';
   s0t(nx+i,:) = (alphahat2+Phat2*Bz2'*inv(Ft2)*nu2')'; 
   s0t(2*nx+i,:) = (alphahat3+Phat3*Bz3'*inv(Ft3)*nu3')'; 
   s0t(3*nx+i,:) = (alphahat4+Phat4*Bz4'*inv(Ft4)*nu4')'; 
   
   ltrvecP0(i,:)=(ltrvec(Phat1-Phat1*Bz1'*inv(Ft1)*Bz1*Phat1'))'; 
   ltrvecP0(nx+i,:) = (ltrvec(Phat2-Phat2*Bz2'*inv(Ft2)*Bz2*Phat2'))'; 
   ltrvecP0(2*nx+i,:) = (ltrvec(Phat3-Phat3*Bz3'*inv(Ft3)*Bz3*Phat3'))'; 
   ltrvecP0(3*nx+i,:) = (ltrvec(Phat4-Phat4*Bz4'*inv(Ft4)*Bz4*Phat4'))'; 
    
   i=i+1;
   
 end 
   
   Ft = Ft -yhat*yhat'; 
   nuhat(t,:) = data(t,:)-yhat';
   
   % llik_c computation 
   
   llik_htmax = max(llik_ht); 
   llik_c(t) = llik_c(t) + llik_htmax + log(probx2'*exp(llik_ht-llik_htmax)); 
   
   % updating mixture probabilities and the state probabilities 
   
   probx2 = (probx2.*exp(llik_ht-llik_htmax))/(probx2'*exp(llik_ht-llik_htmax));
   probs  = [sum(probx2(1:nx));sum(probx2(nx+1:2*nx));sum(probx2(2*nx+1:3*nx));sum(probx2(3*nx+1:4*nx))]; 
   
   % collapse the state vector : Kim and Nelson (1999)'s approximation 
   
   zeroprob = (probx2<=1e-5); 
   
   if sum(zeroprob) >= (4*nx-nxmax)
   % simply delete the zero probability states 
   s1t = s0t(zeroprob==0,:);
   ltrvecP1 = ltrvecP0(zeroprob==0,:); 
   probx = probx2(zeroprob==0); 
   
   else

       % collapse some of the states
   selmat       = kron(eye(nx),[1 1 1 1]); 
   sumzeroprob  = selmat*zeroprob;
   del          = sumzeroprob==4;  
   nx4          = sum(del); 
   
   s1t = zeros(nx-nx4,ns);
   ltrvecP1 = zeros(nx-nx4,ntr); 
   probx = zeros(nx-nx4,1); 
   j=1; 
   
   i=1;
   while i<=nx; 
   
       i2 = (i-1)*4+1; 
       
       if sumzeroprob(i)<4 ; 
          % eliminate histories if all the subhistories have insignificant prob.    
          % collapse density 
            probx(j) = probx2(i2)+probx2(i2+1)+probx2(i2+2)+probx2(i2+3); 
            s1t(j,:) = (probx2(i2)/probx(j))*s0t(i2,:)+(probx2(i2+1)/probx(j))*s0t(i2+1,:)+(probx2(i2+2)/probx(j))*s0t(i2+2,:)+(probx2(i2+3)/probx(j))*s0t(i2+3,:);
            ltrvecP1(j,:)=ltrvec((probx2(i2))/(probx(j))*(ltrmat(ltrvecP0(i2,:)')+s0t(i2,:)'*s0t(i2,:)) +...
                          (probx2(i2+1))/(probx(j))*(ltrmat(ltrvecP0(i2+1,:)')+s0t(i2+1,:)'*s0t(i2+1,:))+...
                          (probx2(i2+2))/(probx(j))*(ltrmat(ltrvecP0(i2+2,:)')+s0t(i2+2,:)'*s0t(i2+2,:))+... 
                          (probx2(i2+3))/(probx(j))*(ltrmat(ltrvecP0(i2+3,:)')+s0t(i2+3,:)'*s0t(i2+3,:))+...
                           -s1t(j,:)'*s1t(j,:))';
            j=j+1;
        
       end
       i=i+1;
   end
   end
   probx = probx./sum(probx);
   nnx = size(s1t,1);
      if nnx>nxmax;
       disp('error : nx is greater than nxmax')
       t=nobs+1; 
      end
     % store information for future use 
   probsT(t,:) = probs'; 
t=t+1;
end

logl = sum(llik_c);
elseif eu>0  
logl = -1000000;
 end                  % determinacy case 
end                   %bound 

    
% Backward Simulation 
t = nobs;
ssT_draw(t,:) = probsT(t,:);

while t>1 
 t = t-1; 
 probs1 = ssT_draw(t+1,:)*trans(:,1)*probsT(t,1)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 probs2 = ssT_draw(t+1,:)*trans(:,2)*probsT(t,2)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 probs3 = ssT_draw(t+1,:)*trans(:,3)*probsT(t,3)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 
 
 ssT_draw(t,:) = [probs1, probs2, probs3, 1-probs1-probs2-probs3]; 
end   
  
 elseif mprior == 5

% parameter loading 

alpha1 = para(1);
alpha2 = para(2);
gamma1 = para(3);
gamma2 = para(4);
beta1  = para(5);
kappa1 = para(6);
tau1   = para(7);
rho_a = para(8);
rho_f = para(9);
rho_i  = para(10);
sd_a1  = para(11);
sd_f1  = para(12);
sd_i1   = para(13);
sd_a2  = para(14);
sd_f2  = para(15);
sd_i2   = para(16);
gy     = para(17);
lnA_0  = para(18);
yst    = para(19);
p11 = para(20);
p22 = para(21);
p11_a = para(22);
p22_a = para(23);
lnP_0 = para(24);
sd_p1 = para(25);
sd_p2=sd_p1;

bound  = zeros(length(para),2);
bound(18:19,1) = -100; 
bound(:,2) = [100;100;100;100;0.9999;100;100;0.9999;0.9999;0.9999;100;100;100;100;100;100;100;100;100;1;1;1;1;100;100]; 
     
crit1 = para-bound(:,1); 
crit2 = bound(:,2)-para; 

        if min(crit1)<0 || min(crit2)<0 
logl = -10000000; 
           return
        else

P_m  = [p11 1-p11;
       1-p22 p22];
   
P_a = [p11_a 1-p11_a;
       1-p22_a p22_a];   

P = kron(P_m ,P_a);



iss    = gy-log(beta1);

%----------------------------------------------------------------------
% check the roots of the system using gensys
%----------------------------------------------------------------------
A = [p11 1-p11 (p11/tau1) ((1-p11)/tau1);
     1-p22 p22 ((1-p22)/tau1) (p22/tau1);
     0 0 beta1*p11 beta1*(1-p11);
     0 0 beta1*(1-p22) beta1*p22];

B = [1+(gamma1/tau1) 0 (alpha1/tau1) 0;
     0 1+(gamma2/tau1) 0 (alpha2/tau1);
     -kappa1 0 1 0;
     0 -kappa1 0 1];

eu =0; 
if (min(abs(real(eig(inv(A)*B))))) < 1
    disp('*************  DL determinacy condition violated  *************')
    eu =1; 
end

if eu ==0 
%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------
A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
            -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];           
solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
            -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];         
solvec_f = inv(A_Af_mat) * [0;0;kappa1;kappa1;0;0];

A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
            -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
            -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
            0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;                      
            -gamma1, 0, -alpha1, 0, 1, 0;
            0, -gamma2, 0, -alpha2, 0, 1];        
solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];
 y1a = solvec_a(1);
 y2a = solvec_a(2);
 pi1a = solvec_a(3);
 pi2a = solvec_a(4);
 i1a = solvec_a(5);
 i2a = solvec_a(6);
 y1f = solvec_f(1);
 y2f = solvec_f(2);
 pi1f = solvec_f(3);
 pi2f = solvec_f(4);
 i1f = solvec_f(5);
 i2f = solvec_f(6);
 y1i = solvec_i(1);
 y2i = solvec_i(2);
 pi1i = solvec_i(3);
 pi2i = solvec_i(4);
 i1i=solvec_i(5);
 i2i=solvec_i(6);

%----------------------------------------------------------------------
% compute coefficients for measurement equations 
%----------------------------------------------------------------------

rho1     = diag([rho_a;rho_f;rho_i;1;1]);
rho1(4,1)= rho_a; 


TC  = [0;0;0;gy;0]; 

sd1 = [sd_a1;sd_f1;sd_i1]; 
sd1  = [sd1;sd_a1;sd_p1]; sd2=[sd_a2;sd_f2;sd_i2;sd_a1;sd_p2]; sd3=sd1; sd4=sd2;



% system matrices for the macro part 
mac1 = [y1a y1f y1i 1 0; pi1a pi1f pi1i 0 1;i1a i1f i1i 0 1]; 
mac2 = [y2a y2f y2i 1 0; pi2a pi2f pi2i 0 1;i2a i2f i2i 0 1]; 

% Markov Chain transition probability 
trans = P';  

% compute the stationary distribution 
A = [eye(length(P))-trans;ones(1,length(P))]; 
probs = inv(A'*A)*(A'*[zeros(length(P),1);1]); 
Vs1 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0;0]);
Vs2 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0;0]);
Vs3 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2);0;0]);
Vs4 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2);0;0]);


Vs =probs(1)*Vs1+probs(2)*Vs2+probs(3)*Vs3+probs(4)*Vs4; 


 % Inflation Persistence 
sumw1      = pi1a^2*Vs1(1,1)+pi1f^2*Vs1(2,2)+pi1i^2*Vs1(3,3); 
sumw2      = pi1a^2*Vs2(1,1)+pi1f^2*Vs2(2,2)+pi1i^2*Vs2(3,3);
sumw3      = pi2a^2*Vs3(1,1)+pi2f^2*Vs3(2,2)+pi2i^2*Vs3(3,3); 
sumw4      = pi2a^2*Vs4(1,1)+pi2f^2*Vs4(2,2)+pi2i^2*Vs4(3,3);


weight11   = (pi1a^2)*Vs1(1,1)/sumw1; 
weight12   = (pi1f^2)*Vs1(2,2)/sumw1; 
weight21   = (pi1a^2)*Vs2(1,1)/sumw2; 
weight22   = (pi1f^2)*Vs2(2,2)/sumw2; 
weight31   = (pi2a^2)*Vs3(1,1)/sumw3; 
weight32   = (pi2f^2)*Vs3(2,2)/sumw3; 
weight41   = (pi2a^2)*Vs4(1,1)/sumw4; 
weight42   = (pi2f^2)*Vs4(2,2)/sumw4; 

rho_inf1   = rho_a*weight11+rho_f*weight12+rho_i*(1-weight11-weight12);  
rho_inf2   = rho_a*weight21+rho_f*weight22+rho_i*(1-weight21-weight22); 
rho_inf3   = rho_a*weight31+rho_f*weight32+rho_i*(1-weight31-weight32);  
rho_inf4   = rho_a*weight41+rho_f*weight42+rho_i*(1-weight41-weight42); 

rho_inf    = [rho_inf1;rho_inf2;rho_inf3;rho_inf4;]; 


% auxiliary parameters 
nobs = size(data,1);
ns   = size(rho1,1);   % dimension of state vector
nz   = size(data,2);      % dimension of observed variables 
ntr  = ns*(ns+1)/2;
nuhat = zeros(nobs,nz); 
nhist = 2;              
nxmax = 4^nhist; 

Az      = [ yst ; 0 ; iss]; 
Bz1    = mac1; Bz2 = mac2; 


% The matrix for storing information from the forward simulation 
probsT = zeros(nobs,4); 
ssT_draw = zeros(nobs,4); 
probsT(1,:) = probs';


% condition on the first observation : assume the first regime 
s00      = zeros(ns,1); 
s00(4)   = lnA_0; 
s00(5)   = lnP_0; 

llik_c   = zeros(nobs,1);       % log likelihood contribution
nu        = data(1,:)'-(Az+Bz1*(TC+rho1*s00)); 
Ft        = Bz1*Vs*Bz1'; 
llik_c(1) =  -0.5*nz*log(2*pi)-0.5*log(det(Ft))-0.5*nu'*inv(Ft)*nu;

s00      = s00+(Vs*Bz1')*inv(Ft)*nu; 
ltrvecP1 = (ltrvec(Vs-Vs*Bz1'*inv(Ft)*(Vs*Bz1')'))';
probx    = 1;
s1t = s00'; 
t=2;
% log likelihood by Kalman filter
while t<nobs+1
nx    = size(s1t,1);   
probs = trans*probs; 
probx2 = kron(probs,probx); 

llik_ht  = zeros(4*nx,1);      
s0t      = zeros(4*nx,ns);    
ltrvecP0 = zeros(4*nx,ntr);   
yhat     = zeros(nz,1); 
Ft       = zeros(nz,nz); 

i = 1; 
 while i<=nx; 

     
   % forecasting
   
    alphahat = TC+rho1*s1t(i,:)';
   P1t      = ltrmat((ltrvecP1(i,:))');
   Phat1    = rho1*P1t*rho1'+(diag(sd1))^2;
   Phat2    = rho1*P1t*rho1'+(diag(sd2))^2;
   Phat3    = rho1*P1t*rho1'+(diag(sd3))^2; 
   Phat4    = rho1*P1t*rho1'+(diag(sd4))^2; 
  
  
   yhat1    = Az  + Bz1*alphahat;
   yhat2    = Az  + Bz2*alphahat; 
      
   nu1      = data(t,:)-yhat1'; 
   nu2      = data(t,:)-yhat2'; 
    
   Ft1      = Bz1*Phat1*Bz1'; Ft2 = Bz1*Phat2*Bz1'; Ft3 = Bz2*Phat3*Bz2';  Ft4 = Bz2*Phat4*Bz2';  
  

   yhat     = yhat + (probx2(i)+probx2(nx+i))*yhat1+(probx2(2*nx+i)+probx2(3*nx+i))*yhat2;
  
   Ft       = Ft+probx2(i)*(Ft1+yhat1*yhat1')+probx2(nx+i)*(Ft2+yhat1*yhat1')+probx2(2*nx+i)*(Ft3+yhat2*yhat2')+probx2(3*nx+1)*(Ft4+yhat2*yhat2');

   % compute likelihood 
   
   llik_ht(i,1)    = -0.5*nz*log(2*pi)-0.5*log(det(Ft1))-0.5*nu1*inv(Ft1)*nu1'; 
   llik_ht(nx+i,1)   = -0.5*nz*log(2*pi)-0.5*log(det(Ft2))-0.5*nu1*inv(Ft2)*nu1'; 
   llik_ht(2*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft3))-0.5*nu1*inv(Ft3)*nu2'; 
   llik_ht(3*nx+i,1) = -0.5*nz*log(2*pi)-0.5*log(det(Ft4))-0.5*nu1*inv(Ft4)*nu2'; 
  

 
   % updating step 
   s0t(i,:)    = (alphahat+Phat1*Bz1'*inv(Ft1)*nu1')'; s0t(nx+i,:) = (alphahat+Phat2*Bz1'*inv(Ft2)*nu1')'; 
   s0t(2*nx+i,:)= (alphahat+Phat3*Bz2'*inv(Ft3)*nu2');  s0t(3*nx+i,:) = (alphahat+Phat4*Bz2'*inv(Ft4)*nu2'); 
  

   ltrvecP0(i,:)=(ltrvec(Phat1-Phat1*Bz1'*inv(Ft1)*Bz1*Phat1'))';     ltrvecP0(nx+i,:)   = (ltrvec(Phat2-Phat2*Bz1'*inv(Ft2)*Bz1*Phat2'))'; 
   ltrvecP0(2*nx+i,:)=(ltrvec(Phat3-Phat3*Bz2'*inv(Ft3)*Bz2*Phat3'))';  ltrvecP0(3*nx+i,:) = (ltrvec(Phat4-Phat4*Bz2'*inv(Ft4)*Bz2*Phat4'))';   
  

   i=i+1;
     
 end 
   
   Ft = Ft -yhat*yhat'; 
   nuhat(t,:) = data(t,:)-yhat';
   
   % llik_c computation 
   
   llik_htmax = max(llik_ht); 
   llik_c(t) = llik_c(t) + llik_htmax + log(probx2'*exp(llik_ht-llik_htmax)); 
   
   % updating mixture probabilities and the state probabilities 
   
  probx2 = (probx2.*exp(llik_ht-llik_htmax))/(probx2'*exp(llik_ht-llik_htmax));
   probs  = [sum(probx2(1:nx));sum(probx2(nx+1:2*nx));sum(probx2(2*nx+1:3*nx));sum(probx2(3*nx+1:4*nx))];   
  


   
   % collapse the state vector : Kim and Nelson (1999)'s approximation 
   
   zeroprob = (probx2<=1e-5); 
   
   if sum(zeroprob) >= (16*nx-nxmax)
   % simply delete the zero probability states 
   s1t = s0t(zeroprob==0,:);
   ltrvecP1 = ltrvecP0(zeroprob==0,:); 
   probx = probx2(zeroprob==0); 
   
   else

       % collapse some of the states
   selmat       = kron(eye(nx),[1 1 1 1 ]); 
   sumzeroprob = selmat*zeroprob;    
   del          = sumzeroprob==4;  
   nx16          = sum(del); 
     
   s1t = zeros(nx-nx16,ns);
   ltrvecP1 = zeros(nx-nx16,ntr); 
   probx = zeros(nx-nx16,1); 
   j=1; 
   
   i=1;
   while i<=nx; 
   
       i2 = (i-1)*4+1; 
       
      if sumzeroprob(i)<4 ; 
          % eliminate histories if all the subhistories have insignificant prob.    
          % collapse density 
            probx(j) = sum(probx2(i2:i2+3)); 
  s1t(j,:) = (probx2(i2)/probx(j))*s0t(i2,:)+(probx2(i2+1)/probx(j))*s0t(i2+1,:)+(probx2(i2+2)/probx(j))*s0t(i2+2,:)+(probx2(i2+3)/probx(j))*s0t(i2+3,:);

           ltrvecP1(j,:)=ltrvec((probx2(i2))/(probx(j))*(ltrmat(ltrvecP0(i2,:)')+s0t(i2,:)'*s0t(i2,:)) +...
                          (probx2(i2+1))/(probx(j))*(ltrmat(ltrvecP0(i2+1,:)')+s0t(i2+1,:)'*s0t(i2+1,:))+...
                          (probx2(i2+2))/(probx(j))*(ltrmat(ltrvecP0(i2+2,:)')+s0t(i2+2,:)'*s0t(i2+2,:))+... 
                          (probx2(i2+3))/(probx(j))*(ltrmat(ltrvecP0(i2+3,:)')+s0t(i2+3,:)'*s0t(i2+3,:))-s1t(j,:)'*s1t(j,:))';
    
            j=j+1;
        
       end
      

       i=i+1;
   end
   end
   probx = probx./sum(probx);
   nnx = size(s1t,1);
      if nnx>nxmax;
       disp('error : nx is greater than nxmax')
       t=nobs+1; 
      end
     % store information for future use 
   probsT(t,:) = probs'; 
t=t+1;
end

logl = sum(llik_c);
elseif eu>0  
logl = -1000000;
 end                  % determinacy case 
end                   %bound 

    
% Backward Simulation 
t = nobs;
ssT_draw(t,:) = probsT(t,:);

while t>1 
 t = t-1; 
 probs1 = ssT_draw(t+1,:)*trans(:,1)*probsT(t,1)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 probs2 = ssT_draw(t+1,:)*trans(:,2)*probsT(t,2)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 probs3 = ssT_draw(t+1,:)*trans(:,3)*probsT(t,3)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 probs4 = ssT_draw(t+1,:)*trans(:,4)*probsT(t,4)/(ssT_draw(t+1,:)*trans*probsT(t,:)');
 
 
 ssT_draw(t,:) = [probs1, probs2, probs3,probs4]; 
end   
  
  
end  
